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Introduction 

The  scientific  objective  of  this  research  is  to  develop  and  characterize  an  MRI  coupled  fluorescence  molecular 
tomography  (FMT)  system  and  associated  reconstruction  algorithms  to  image  exogenous  fluorescent  agents  in 
human  breast.  Imaging  the  spatial  distribution  of  unspecific  and/or  targeted  fluorescent  probes  may  reveal 
functional  information  useful  in  determining  whether  a  region  of  interest  is  indeed  malignant,  thus  reducing  the 
high  number  of  false  positives  and  associated  biopsies,  and  aiding  in  treatment  planning.  Fluorescence  imaging 
at  depth  in  tissue  is  a  challenging  problem  and  FMT  images  often  suffer  from  relatively  poor  resolution,  depth 
dependent  sensitivity,  and  inaccurate  quantification,  especially  in  complex  tissue  volumes  such  as  human 
breast.  However,  incorporating  tissue  structural  information  obtained  from  a  high  resolution  clinical  imaging 
modality,  such  as  MRI,  into  the  FMT  image  reconstruction  routine  may  improve  the  spatial  resolution  and 
imaging  accuracy  of  fluorescence  yield  images.  An  FMT  system  incorporated  directly  into  the  MR  magnet 
bore  for  simultaneous  optical  and  MR  image  acquisition  provides  a  unique  opportunity  to  couple  the  high- 
resolution  information  from  MRI  with  the  highly  sensitive  optical  data  to  obtain  accurate  images  of 
fluorescence  drug  distribution  in  complex  domains.  This  system  and  the  relevant  algorithms  have  been 
developed  in  the  funding  period  of  this  grant  and  preliminary  results  indicate  a  dramatic  improvement  using 
“spatial  priors”  from  MR  images.  This  report  summarizes  the  major  accomplishments  to  date. 

Body 

Significant  progress  has  been  made  toward  the  proposed  aims  of  this  project  including  development  and  testing 
of  the  experimental  system  (Aims  1  and  2),  advancement  and  validation  of  the  reconstruction  algorithms  (Aims 
3  and  4),  and  development  of  reconstruction  algorithms  to  image  exogenous  agents  using  the  absorption  profile 
of  the  drug  (Aim  6).  A  description  of  the  primary  accomplishments  follows. 
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Fig.  1.  Schematic  of  the  experimental  system  which  couples  tomography  into  a  Philips  3T 
MRI.  The  newly  developed  component  of  the  system  is  the  spectroscopy  detection  system 
which  consists  of  16  spectrometers  for  fluorescence  molecular  tomography. 
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Imaging  system 
The  fluorescence 
imaging  system  is 
coupled  into  a 
Philips  3T  MRI  and 
is  shown 

schematically  in  Fig. 
1 ,  which  also  depicts 
a  legacy  PMT 
frequency  domain 
system  for 

recovering 
background  tissue 
optical  properties 
necessary  for  the 
fluorescence  image 
reconstruction.  The 
newly  constructed 
“Spectroscopy 
system”  is 

composed  of  sixteen 
Insight 
spectrometers 
(Acton  Research) 


with  Pixis  cooled  CCD  cameras  which  provide  highly  sensitive  spectrally  resolved  detection  of  the  collected 
light.  Each  spectrometer  is  coupled  to  the  patient  interface  via  a  13  meter  bifurcated  silica  fiber  (a  total  16 
fibers)  custom  designed  for  this  system.  These  custom  designed  bifurcated  fibers  circumscribe  the  tissue 
volume  being  imaged  and  are 
used  for  both  source  delivery 
and  detection.  One  end  of  the 
bifurcated  fibers  is  coupled  into 
the  spectrometers  with  custom 
designed  and  manufactured 
input  optics  (Fig.  2)  each  of 
which  contains  an  automated, 
six-position  filter  wheel  to  help 
filter  the  excitation  laser  source 
from  the  fluorescence  signal. 

The  input  optics  collimate  the 
light  from  the  fiber,  pass  it 
through  a  filter  of  the  user’s 
choice,  then  focus  the  light  on 
the  spectrometer  entrance  slit 
matching  the  f-number  of  the 
spectrometer  system. 

Photographs  of  the  experimental 
system  coupled  into  the  MR1 
magnet  bore  are  shown  in  Fig.  3. 

All  spectrometers  are  controlled 

Fig.  2.  Custom  designed  input  optics  match  the  f-number  of  the  input  light  to  that  of 
simultaneously  with  programs  th=  spectrograph 
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Fig.  3.  Photographs  of  the  spectroscopy  system  (left)  in  the  MR  control  room.  The  fibers  run  through  Ports  in  the  wall  into  the  MR  bore  (right). 


written  in  Labview  which  cycle  the  laser  source  sequentially  through  each  fiber  while  the  other  15  fibers  detect 
the  emitted  light  from  the  tissue  volume,  resulting  in  a  total  of  240  data  points  for  a  given  acquisition.  Exposure 
times  are  determined  automatically  during  acquisition  to  optimize  the  SNR.  Preliminary  results  using  this 
system  are  presented  later  in  this  report. 
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Lutetium  Texaphyrin  and  spectral  fitting 

The  normalized  absorption  and  fluorescence  spectra  for  Lutetium  Texaphyrin  (LuTex)  are  shown  in  Fig.  4. 
Because  of  the  large  absorption  peak  at  735  nm,  this  probe  is  ideal  for  comparing  both  absorption  tomography 
and  fluorescence  tomography  (Aim  6).  Similar  absorption/fluorescence  spectra  were  recorded  for  Gadolinium 
Texaphyrin  which  has  a  lower  fluorescence  yield. 


Fig.  4.  The  normalized  absorbance  and  fluorescence 
emission  spectra  of  Lutetium  Texaphyrin  are  shown. 


The  spectroscopic  detection  system  offers  a  unique 
opportunity  to  filter  the  excitation  source  light  from  the 
relatively  weak  fluorescence  emission  signal  beyond  the 
filtering  provided  by  the  interference  filters  in  the  input 
optics  and  the  spectrograph  grating.  A  simple  least 
squares  fitting  algorithm  was  developed  to  fit  the  data  to 
pre-recorded  basis  spectra  of  the  fluorescence  emission 
peak  and  excitation  leakage,  thus  decoupling  the 
fluorescence  and  excitation  signals.  Figure  5  provides  an 
example  of  the  fitting  procedure.  The  algorithm 
determines  the  amount  of  each  of  the  pre-recorded  basis 
spectra  depicted  in  (b)  used  to  fit  to  the  data  in  (a).  The 
integrated  fluorescence  basis  function  is  then  used  as  the 
fluorescence  intensity  in  the  reconstruction  algorithm. 
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CCD  Hid  lumber 
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CCD  Hid  lumber 


Fig.  5.  An  example  of  a  pair  of  basis  spectra  for  the  excitation  and  fluorescence 
light  (in  counts/s  as  a  function  of  CCD  pixel  number)  is  shown  in  (b).  These 
spectra  are  recorded  for  each  detector  prior  to  imaging.  In  practice,  the  basis 
spectra  are  used  to  perform  a  least  squares  fit  (a)  to  the  spectrum  measured  for 
each  source- detector  pair  to  determine  the  relative  contribution  of  the 
fluorescence  and  excitation  light  to  the  measured  response. 


Spatial  priors  image  reconstruction  implementation 

The  relatively  poor  spatial  resolution  characteristic  of  diffuse  optical  tomography  (DOT)  has  previously  been 
addressed  by  incorporating  spatial  prior  information  from  other  imaging  modalities[l-4]  but  has  not  been 
previously  applied  to  fluorescence  tomography.  The  approach  presented  here  incorporates  spatial  information 
from  MR  images  into  the  reconstruction  algorithm  using  a  Laplacian-type  matrix  to  regularize  the  inversion 
formulation  as  described  in  Ref. [5],  Results  of  this  approach  with  simulated  and  phantom  data  indicate 
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dramatic  improvements  in  imaging  performance,  especially  for  complex  tissue  volumes.  Figure  6  depicts  a  T1 
weighted  MR  image  of  a  human  breast  acquired  in  previous  studies  with  our  frequency  domain  diffuse  optical 
tomographic  (DOT)  system.  This  image  was  segmented  into  adipose  and  fibro-glandular  regions  and  a  third 
region  was  added  to  simulate  a  malignant  tumor.  Each  region  was  tagged  with  clinically  relevant 
concentrations  of  endogenous  chromophores  (oxy  and  deoxy  hemoglobin,  water)  as  well  as  estimated 
concentrations  for  Lutetium  Texaphyrin,  one  of  the  drugs  currently  under  investigation.  The  domain  was 
meshed  and  simulated  noisy  data  of  excitation  and  fluorescence  emission  were  generated.  The  data  was  then 
used  to  reconstruct  images  of  the  background  optical  properties  and  fluorescence  yield  using  the  non-spatial 
algorithms  and  the  Laplacian-type  spatial  priors  implementation. 


Fig.  6.  An  axial  T1  weighted  MR  slice  of  a  human  breast  is  shown  in 
(a),  from  which  the  test  domain  for  simulation  studies  was  derived. 

Darker  regions  indicate  fibro-glandular  tissue  imbedded  in  adipose 
tissue  indicated  by  lighter  values.  The  image  was  acquired  during  a 
clinical  exam  with  one  of  our  MR-coupled  NIR  tomography  systems 
and  shows  the  indentations  caused  by  the  fiber  optic  probes.  The 
test  domain  (b)  was  a  discretization  of  the  MR  image  thresholded  into 
regions.  The  yellow  anomaly  was  added  to  simulate  a  targeted 
cancerous  tumor. 

Recovered  images  are  shown  in  Fig.  7  for  a  tumor  located  near  the  edge  (left)  and  near  the  center  (right)  of  the 
breast  volume.  The  bottom  row  of  images  show  the  values  of  fluorescence  yield  and  demonstrate  no  ability  to 
quantify  the  tumor  region  without  using  spatial  priors.  A  dramatic  improvement  is  revealed  using  the  spatial- 
priors  approach  which  provides  qualitatively  and  quantitatively  accurate  images  of  fluorescence  yield. 
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Fig.  7.  Target  and  recovered  values  for  reconstruction  implementations  using  no  prior  information  and  with  spatial  prior 
information. 


The  spatial  priors  implementation  was  also  validated  with  experimental  data  acquired  using  the  system 
described  above.  A  55mm  diameter  solid  phantom  with  a  drilled  well  was  used  for  these  experiments  and  is 

depicted  in  Fig.  8.  Different  concentrations  of  the  fluorophore  Lutetium 
Texaphyrin  in  intralipid  and  water  acted  as  the  fluorescent  inclusion  in  this 
simple  test  domain.  A  690  nm  CW  laser  source  was  used  to  excite  the 
fluorophore  and  the  multi-channel  spectrometer  system  used  to  detect  the 
emitted  fluorescence.  Image  acquisition  time  totalled  4  minutes  for  these 
experiments  and  further  details  are  presented  in  Ref. [5],  Images  were 
reconstructed  with  both  the  non-spatial  priors  and  spatial  priors 
implementations  and  the  results  presented  in  Fig.  9.  The  spatial  priors 
implementation  produces  more  accurate  reproductions  of  the  fluorescent 
inclusion  and  suppresses  edge  artifacts  and  image  noise.  Further 
characterization  of  the  experimental  system  including  imaging  larger,  more 
complicated  volumes  will  be  studied  as  part  of  the  second  year  of  funding. 


Fig.  8.  Solid  phantom  used  for 
experimental  validation. 
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Fig.  9.  Recovered  images  of  fluorescence  yield  are  shown  for  varying  concentrations  of  Lutex.  The  14  mm  diameter  fluorescent 
inclusion  was  embedded  in  a  55  mm  diameter  solid  epoxy  tissue  simulating  phantom.  Images  were  generated  from  the  same  data  using 
algorithms  based  on  no-priors  and  spatial  soft  prior  implementations. 


Absorption  tomography  for  imaging  exogenous  contrast 

Another  approach  to  image  exogenous  contrast  uses  the  optical  absorption  characteristics  of  a  contrast  agent 
rather  than  the  emitted  light  from  a  fluorescent  contrast  agent.  One  of  the  stated  aims  of  this  funded  research  is 
to  compare  the  absorption  and  fluorescence  based  approaches  to  image  the  spatial  distribution  of  an  exogenous 
drug.  A  “spectrally  constrained”  approach  which  recovers  images  of  endogenous  tissue  chromophore  (oxy  and 
deoxy  hemoglobin,  water)  concentrations  and  scattering  parameters  using  the  chromophores’  known  spectral 
absorption  profiles  as  prior  information  has  been  shown  in  previous  studies  to  accurately  image  the  endogenous 
chromophore  concentration  in  human  breast[l,  2,  6].  Under  the  current  grant,  this  approach  has  been  expanded 
to  incorporate  exogenous  absorbing  drugs.  Simulation  studies  have  been  completed  which  demonstrate 
simultaneous  recovery  of  endogenous  and  exogenous  chromophores  and  scattering  parameters  using  multi¬ 
wavelength  spectrally-constrained  absorption  tomography.  The  target  values  of  endogenous  and  exogenous 
Lutex  chromophore  concentrations  depicted  in  Fig.  10  were  used  to  generate  simulated  noisy  frequency  domain 
data  at  ten  discrete  wavelengths.  This  data  was  used  to  complete  a  spectrally  constrained  reconstruction  to 
recover  images  of  the  target  domain.  The  results  presented  in  Fig.  10  indicate  reasonable  qualitative  and 
quantitative  accuracy  for  the  recovery  of  both  the  endogenous  chromophores  as  well  as  LuTex.  This  approach 
will  be  formally  compared  to  the  fluorescence  imaging  case  in  the  second  year  of  the  funding  period. 
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Fig.  10.  Recovered  and  target  images  of  endogenous  and  exogenous  chromophores  in  a  realistic  breast  volume.  Images  were 
reconstructed  using  a  spectrally  constrained  technique  which  incorporates  the  spectral  absorption  features  as  prior  knowledge. 


Key  research  accomplishments 

•  Constructed  experimental  system  and  conducted  initial  calibration  studies. 

•  Acquired  preliminary  data  and  verified  the  feasibility  of  imaging  fluorescence  yield  down  to  300nM  of 
Lutex. 

•  Incorporated  tissue  structural  information  obtained  from  MR  images  into  the  image  reconstruction 
algorithms  and  demonstrated  a  dramatic  qualitative  and  quantitative  improvement  in  imaging 
performance. 

•  Developed  and  tested  algorithms  for  simultaneous  reconstruction  of  endogenous  and  exogenous 
chromophores. 

Reportable  outcomes 

The  first  year  of  this  traineeship  grant  has  led  to  one  first  author  publication  manuscript  currently  under  review 
and  three  first  author  conference  presentations. 

Scott  C.  Davis,  Hamid  Dehghani,  Jia  Wang,  Shudong  Jiang,  Brian  W.  Pogue  and  Keith  D.  Paulsen,  “Image- 
guided  diffuse  optical  fluorescence  tomography  implemented  with  Laplacian-type  regularization”,  Optics 
Express,  (in  review). 

Scott  C.  Davis,  Hamid  Dehghani,  Jia  Wang,  Shudong  Jiang,  Brian  W.  Pogue  and  Keith  D.  Paulsen,  “MR 
guided  reconstruction  of  Fluorescence  Diffuse  Optical  Tomography”,  Oral  presentation  at  the  International 
Society  for  Optical  Engineering  (SPIE)  BiOS-2007  in  Photonics  West-2007,  San  Jose,  California,  January  20- 
25,  2007. 

Scott  C.  Davis,  Shudong  Jiang,  Hamid  Dehghani,  Roger  Springett,  Keith  D.  Paulsen,  Brian  W.  Pogue,  “Multi- 
spectral,  MR-guided  Diffuse  Optical  Tomography  for  Imaging  Fluorescence  and  Absorption  Exogenous 
Contrast”,  Poster  presentation  at  the  Fifth  Annual  Meeting  of  the  Society  for  Molecular  Imaging,  Hilton 
Waikoloa  Village,  Hawaii,  USA,  August  30  -  September  2,  2006. 
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Scott  C.  Davis,  Brian  W.  Pogue,  Subhadra  Srinivasan,  Hamid  Dehghani  and  Keith  D.  Paulsen,  “Development 
of  spectrally-constarined  diffuse  optical  tomography  for  imaging  exogenous  contrast  agents”,  Biomedical 
Optics  Topical  Meeting  on  CD-ROM,  The  Optical  Society  of  America,  Washington  DC,  March  19-22,  2006. 
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Conclusions 

The  funding  period  covered  by  this  report  has  produced  a  working,  validated  experimental  fluorescence 
tomography  system  that  uses  a  parallel  spectrometer-base  detection  system  which  couples  into  a  clinical  MRI 
magnet.  Preliminary  results  using  this  system  demonstrate  qualitatively  accurate  images  down  to  300nm  of 
Lutex  in  a  simple  solid/liquid  tissue  phantom  using  the  newly  developed  spatial-priors  reconstruction 
implementation.  In  simulation  studies  of  realistic  breast  domains,  spatial  priors  provide  dramatic  improvements 
in  image  quality  and  quantification.  Several  conference  presentations  and  a  journal  manuscript  have  been 
produced  from  these  encouraging  results.  Additionally,  a  spectrally-constrained  algorithm  which  recovers  the 
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concentration  of  exogenous  absorption  contrast  agents  was  developed  and  results  from  simulated  realistic  breast 
volumes  indicate  a  promising  ability  to  recover  contrast  agent  concentration. 
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Abstract:  A  promising  method  to  incorporate  tissue  structural  information 
into  the  reconstruction  of  diffusion-based  fluorescence  imaging  is 
introduced.  The  method  regularizes  the  inversion  problem  with  a 
Laplacian-type  matrix,  which  inherently  smoothes  pre-defmed  tissue,  but 
allows  discontinuities  between  adjacent  regions.  The  technique  is  most 
appropriately  used  when  fluorescence  tomography  is  combined  with 
structural  imaging  systems.  Phantom  and  simulation  studies  were  used  to 
illustrate  significant  improvements  in  quantitative  imaging  and  linearity  of 
response  with  the  new  algorithm.  Images  of  an  inclusion  containing  the 
fluorophore  Lutetium  Texaphyrin  (Lutex)  embedded  in  a  cylindrical 
phantom  are  more  accurate  than  in  situations  where  no  structural 
information  is  available,  and  edge  artifacts  which  are  normally  prevalent 
were  almost  entirely  suppressed.  Most  importantly,  spatial  priors  provided 
a  higher  degree  of  sensitivity  and  accuracy  to  fluorophore  concentration, 
though  both  techniques  suffer  from  image  bias  caused  by  excitation  signal 
leakage.  The  use  of  spatial  priors  becomes  essential  for  accurate  recovery  of 
fluorophore  distributions  in  complex  tissue  volumes.  Simulation  studies 
revealed  an  inability  of  the  “no-priors”  imaging  algorithm  to  recover  Lutex 
fluorescence  yield  in  domains  derived  from  T1  weighted  images  of  a  human 
breast.  The  same  domains  were  reconstructed  accurately  to  within  75%  of 
the  true  values  using  prior  knowledge  of  the  internal  tissue  structure.  This 
algorithmic  approach  will  be  implemented  in  an  MR-coupled  fluorescence 
spectroscopic  tomography  system,  using  the  MR  images  for  the  structural 
template  and  the  fluorescence  data  for  region  quantification. 

©  2007  Optical  Society  of  America 
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1.  Introduction 

Imaging  the  spatial  distribution  of  fluorescence  activity  at  depth  in  tissue  is  a  challenging 
problem  that  has  yet  to  impact  clinical  practice.  The  most  popular  approach  is  diffuse  optical 
fluorescence  tomography  (DOFT),  a  model-based  method  which  approximates  photon 
propagation  as  a  diffuse  field  and  computationally  matches  model  parameters  to  measured 
boundary  data.  The  theoretical  framework  is  well  established  [1-5]  and  feasibility  studies 
have  demonstrated  fluorescence  yield  imaging  in  various  phantom  geometries  [6-9].  Much  of 
the  most  recent  research  has  focused  on  high-throughput,  full  body  imaging  of  small  animals 
[10,  11]  but  no  human  images  have  been  published  to  date.  DOFT  produces  low  resolution 
images  compared  to  standard  clinical  modalities  such  as  x-ray  CT  and  MRI  due  to  the  highly 
scattered  photon  fields  and  relatively  sparse  measurement  sampling  of  the  tissue  volume.  The 
reconstruction  problem  is  ill-posed  and  underdetermined  and  large  heterogeneity  in  tissue 
optical  properties  caused  by  complex  tissue  morphology  further  challenges  the  imaging 
algorithm.  Sensitivity  drops  with  increasing  depth  resulting  in  a  non-linear  responsivity 
across  the  imaging  field,  making  it  more  difficult  to  image  larger  tissue  volumes. 
Experiments  in  larger  volumes  that  more  closely  resemble  a  human  breast  typically  consider 
unrealistically  high  fluorophore  contrasts  and  simple  geometries  [9,  12].  Though  these  studies 
are  important  for  advancing  the  understanding  of  the  modality,  improving  image  resolution 
and  contrast  sensitivity  is  critical  for  identifying  a  clinical  role  for  DOFT. 

DOFT  is  an  extension  of  the  more  widely  studied  diffuse  optical  tomography  (DOT) 
which  suffers  similarly  from  low  resolution  and  depth-dependent  contrast  sensitivity  in  the 
absence  of  spatial  or  spectral  prior  information  to  guide  the  solution.  However,  methods  to 
incorporate  highly-resolved  anatomical  data  obtained  from  standard  clinical  modalities  have 
improved  the  quantification  of  the  optically  derived  images  [13-15].  These  hybrid  approaches 
lead  to  a  conceptually  new  application  of  optical  tomography,  one  in  which  the  highly 
resolved  imaging  system  provides  a  structural  template  upon  which  volumetric  optical 
spectroscopic  images  are  constructed.  This  framework  may  be  applied  to  either  absorption 
and  scatter  spectroscopy,  or  fluorescence  spectroscopy.  Additional  challenges  specific  to  the 
fluorescence  case  include  lower  signal  intensity,  excitation  source  contamination  of  the 
fluorescence  emission  measurements  due  to  filter  leakage,  and  tissue  optical  property  effects 
on  both  the  excitation  and  emission  photon  propagation.  Applying  spatial  guidance 
techniques  to  the  more  complicated  DOFT  problem  may  yield  even  larger  gains  in  imaging 
capability. 

This  paper  introduces  a  method  to  dramatically  improve  fluorescence  imaging  at  depth  in 
tissue  by  coupling  MRI  and  fluorescence  tomography.  Tissue  structural  information 
determined  from  standard  T1  and  T2  MR  images  is  encoded  as  a  spatial  filter  in  the  DOFT 
reconstruction  algorithm  and  used  to  guide  the  recovery  of  fluorescence  activity.  In  this 
implementation,  reconstruction  parameters  are  loosely  grouped  into  regions  based  on  tissue- 
type  determined  from  the  MR  images  but  are  permitted  to  update  independently,  giving  rise  to 
the  term  “soft”  spatial  priors  [16].  The  algorithm  is  tested  with  phantom  data  of  Lutetium 
Texaphyrin,  a  near- infrared  photosensitizer  shown  to  accumulate  preferential  in  a  variety  of 
malignancies  [17-20],  recorded  from  a  newly  developed  MR-coupled  spectroscopy-based 
fluorescence  tomography  system.  A  simulation  study  based  on  realistic  geometries  generated 
from  MR  images  of  a  normal  human  breast  serves  as  an  initial  illustration  of  expected 
improvements  provided  by  the  spatial  priors  approach  in  complex  domains. 


2.  Theory 

Assuming  highly  diffuse  media,  the  transport  of  light  at  a  given  modulation  frequency  co  in 
the  presence  of  fluorescence  generated  by  an  external  source  at  the  excitation  wavelength 
of  the  fluorescing  agent,  is  modeled  by  two  diffusion  equations,  where  the  solution  of  the  first 
equation  provides  the  driving  source  term  of  the  second  [5]: 

-V-Kx(r)V®x(r,(o)  +  ({ia  (>')  +  -^-)<&x(r, co)  =  q0(r, co)  (1) 

c(r) 

-  V  •  Km (r)V<Dm(r, co)  +  (ji  (r)  +  -^-)® m(r, co)  =  O x(r, co)W  (r) (2) 
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where  subscripts  x  and  m  represent  the  excitation  and  emission  fluence  at  wavelengths  Xx  and 

Xm,  respectively.  The  intrinsic  optical  parameters  jUa  and  jUs  are  the  absorption  and 

reduced  scattering  coefficients  respectively,  q0(r,Cd)  is  an  isotropic  source  and 
®xm  (/%&>)  is  the  photon  fluence  rate  at  position  r.  The  diffusion  coefficient  is  given  by 

1 

K x  m  = - ; -  and  c(r)  is  the  speed  of  light  in  the  medium  at  any  point,  defined 

3< ) 

by  Cf/n(r),  where  n(r)  is  the  index  of  refraction  at  the  same  location  and  cQ  is  the  speed  of  light 
in  vacuum.  The  fluorescence  parameters  are  the  lifetime  r(r)  and  the  fluorescence  yield 

/7//a  (r),  the  latter  a  product  of  the  fluorophore’s  quantum  efficiency  rj  and  its  absorption 
coefficient  jUa^  (r) . 

The  most  appropriate  description  of  the  air-tissue  boundary  is  derived  with  an  index- 
mismatched  type  III  condition  (also  known  as  Robin  or  mixed),  in  which  some  fraction  of  the 
fluence  at  the  external  boundary  of  the  tissue  exits  and  does  not  return.  The  flux  leaving  the 
external  boundary  is  equal  to  the  fluence  rate  at  the  boundary  weighted  by  a  factor  that 
accounts  for  the  internal  reflection  of  light  back  into  the  tissue.  This  relationship  is  described 
in  the  following  equation: 

<t„,fc®)  +  2^n  (3) 

where  %  is  a  point  on  the  external  boundary,  and  A  depends  upon  the  relative  refractive  index 
(RI)  mismatch  between  tissue  Q  and  air.  A  can  be  derived  from  Fresnel’s  law: 

2/(l-i?n)-l  +  |cos#J3 

A  =  — - t - k — —  (4) 

1  -  |cos  0C  | 

where  6C  =  arcsin (nAIR  /nx )  ,  the  angle  at  which  total  internal  reflection  occurs  for  photons 

(«i  /  nAIR  -  if 

moving  from  region  Q  with  RI  n\  to  air  with  RI  nAiR,  and  =  - - —  .  At  the 

(«,/«„« +i) 

external  boundaries,  nAIR  =  1 . 

2.1  Finite  element  implementation: 

When  the  refractive  index  is  homogenous  (assumed  to  be  1.33  [21]),  the  finite  element 
discretization  of  a  volume  Cl  can  be  obtained  by  subdividing  the  domain  into  D  elements 
joined  at  V  vertex  nodes.  In  the  finite  element  method  (FEM),  m  (r)  is  approximated  by 


the  piecewise  continuous  polynomial  function  ^  Q x mUi(r)Qh  ,  where 

Clh  is  a  finite  dimensional  subspace  spanned  by  basis  functions  {w.(r);z  =  1...F}  chosen 

■W-  h 

to  have  limited  support.  The  problem  of  solving  for  Qxm  becomes  one  of  sparse  matrix 

inversion,  and  in  this  work,  a  bi-conjugate  gradient  stabilized  solver  is  used.  As  developed 
previously  [22,  23],  the  coupled  diffusion  Eqs.  (1)  and  (2)  in  the  FEM  framework  can  be 


expressed  as  a  system  of  linear  algebraic  equations: 

K,  (K)  +  C,  (//„  + F ;  V  =  a  (5) 

l  c(r)  2  A  J 

K.W  +  C.(ji',  +-tt)  +  T7^  K  =-Q.  («) 

ico 

where  the  matrices  Kxm(ic),  Cxm  \fj>a  H - ) ,  and  Fxm  have  entries  given  by: 

c(r) 

(7) 

Q 

cx,m9  =  j  (r)  +  u,  (r>j  (r)d"r  (8) 

V  (9) 

and  the  source  vector  Q0  has  terms 

00, .=\»i(r)9o(r)dnr  (10) 

Q 


The  source  term  is  defined  as  a  Gaussian  distribution,  matching  the  intensity  profile  at  the  tip 
of  the  optical  fiber.  Because  the  source  is  assumed  spherically  isotropic,  modeling  is  more 
accurate  when  it  is  centered  one  scattering  distance  within  the  outer  boundary.  The  source 
vector  Qm  for  fluorescence  re-emission  is  expressed  as 

Qm  =J '«,W  ®x(r, G>)TJMaf (r) 7  \°M%  d"r  (H) 

hi  \  +  [(OT{r)\  \ 

and  is  distributed  throughout  the  domain. 

2.2  The  inverse  model 

2.2.1  No  spatial  priors 

In  the  inverse  problem,  the  goal  is  the  recovery  of  optical  properties  at  each  FEM  node  using 
surface  measurements  of  light  fluence  at  both  the  excitation  and  emission  wavelength 
sequentially  (assuming  the  use  of  two  externally  applied  sources,  one  at  each  wavelength), 
followed  by  fluorescence  yield  reconstruction  at  the  emission  wavelength.  The  computational 

approach  is  to  minimize  the  difference  between  measured  fluence,  ,  at  the  tissue 

surface  and  calculated  data,  from  the  model  Eqs.  (1)  and  (2)  by  adjusting  the  spatial 

distribution  of  the  unknown  parameters  through  minimization  of  the  ‘objective’  function.  The 
objective  function  for  recovering  the  optical  properties  at  the  excitation 
wavelength,  jux  =  (jua  ,  Kx  ) ,  is  given  as 


(12) 


X 


NM  ,  \  NN  (  \ 


i= 1  J=l 

where  NM  is  the  total  number  of  measurements  given  by  the  imaging  system,  NN  is  the 
number  of  parameters  representing  the  optical  property  distribution  which  corresponds  to  the 
number  of  nodes  in  the  reconstruction  mesh,  and  /  is  an  NN  x  NN  identity  matrix.  In  general, 

2  d%2 

%  will  not  equal  zero,  but  the  values  of  jUx  for  which 


is  close  to  zero  can  be 


determined,  based  on  an  initial  estimate  of  JUX  . 
deriving  Newton’s  method,  -  is  evaluated  at 


Following  the  Taylor  series  method  for 
jUx  based  on  an  expansion  around  some 


nearby  point  ju  ,  where  the  second  and  higher  order  terms  are  ignored,  leading  to  the 
iterative  update  equation: 

Ajux  =  [jTJ  +  /l/]"1  JT  (<Df  ‘ ”  -  O'; )  (13) 

where  J  is  the  Jacobian  matrix,  here  calculated  using  the  Adjoint  method  [24].  Equation  (13) 
is  known  as  the  Moore-Penrose  generalized  inverse  and  is  found  to  be  suitable  for  problems 
where  the  number  of  unknowns  to  be  recovered  is  much  larger  than  the  amount  of  information 
(measurements)  available.  In  standard  practice,  /  is  an  identity  matrix,  and  in  this  work  X  is 
some  fixed  fraction  multiplied  by  the  maximum  value  on  the  diagonal  of  the  Hessian  matrix 

J1 J ,  and  is  therefore  updated  at  each  iteration.  To  recover  the  optical  properties  at  the 
emission  wavelength,  jum  =  [jua  ,  K m  ) ,  the  externally  applied  illuminating  source  is  changed 

to  one  at  the  emission  wavelength  and  the  formulation  presented  in  Eq.  (13)  is  used. 

The  recovered  optical  properties  are  used  in  the  fluorescence  yield  reconstruction  which 
also  adheres  to  the  minimization  formulation  presented  in  Eq.  (12).  The  unknown  parameter 

in  Eq.  (13)  becomes  7///a  (r) ,  and  the  Jacobian  can  be  calculated  by  similar  Adjoint 
properties  described  above  and  in  Ref.  [5]. 


2.2.2  Spatial  priors 

Spatial  prior  information  is  incorporated  by  assuming  a  ‘generalized  Tikhonov’  penalty  term 
which  is  similar  in  structure  to  Eq.  (12)  except  that  the  identity  matrix  is  replaced  with  a 
Laplacian-type  matrix,  presented  here  for  the  excitation  field: 

j2=z(®r-4>()  <i4) 

i=\  j=i 

where  NN  is  the  number  of  unknowns  in  the  model  [16,  25,  26].  The  constant,  ) 3  ,  balances 

the  effect  of  the  parameters  with  the  model-data  mismatch  in  the  same  manner  as  X  in  Eq. 
(13).  The  dimensionless  ‘filter’  matrix,  L ,  is  generated  using  MRI-derived  priors  and  its 
construction  is  flexible.  In  this  application,  each  node  in  the  FEM  mesh  is  labeled  according 
to  the  region,  or  tissue  type,  with  which  it  is  associated  (in  the  MR  image).  The  L-matrix 
represents  a  Laplacian-type  structure,  the  diagonal  of  which  is  Lu=l  where  i  is  the  nodal 
index.  When  nodes  i  and  j  are  in  the  same  region  containing  n  nodes,  Luj=-  1/n,  otherwise 
Lij=  0.  This  effectively  relaxes  the  smoothness  constraints  at  the  interface  between  different 
tissues,  in  directions  normal  to  their  common  boundary.  The  effect  on  image  quality  is  similar 
to  that  achieved  through  total  variation  minimization  schemes  [27]  but  easily  encodes  internal 
boundary  information  from  MR  images.  Following  the  same  procedure  as  in  the  no-priors 
case,  the  parameter  update  is  given  by 


^ix=[jtj  +  Pltl]'  jT(<x>Mr  -®,c)  (15) 

where  much  like  X  in  Eq.  (12 ),  j3  is  a  fixed  fraction  multiplied  by  the  maximum  value  on  the 

diagonal  of  JT  J  .  This  update  formulation  is  also  used  to  recover  the  optical  properties  at  the 
emission  wavelength  as  well  as  the  fluorescence  yield,  as  described  in  Section  2.2.1. 

2.2.3  Reconstruction  basis 

A  number  of  different  strategies  for  defining  the  reconstruction  basis  are  possible,  including 
second  mesh  and  pixel  basis  [28,  29].  The  choice  of  reconstruction  basis  allows  for 
computational  efficiency  which  serves  to  reduce  the  number  of  unknowns  in  the  algorithm. 
The  problem  at  hand  is  twofold:  The  forward  problem  requires  that  the  volume  of  interest  is 
subdivided  into  adequate  number  of  sub-domains  which  allow  for  an  accurate  description  of 
the  calculated  fields,  whereas  a  reduction  in  the  number  of  unknowns  improves  the  ill- 
posedness  of  the  problem  in  the  reconstruction  algorithm.  This  is  addressed  by  defining  a 
separate  reconstruction  basis  (different  from  the  meshes  used  in  the  FEM  implementation), 
upon  which  the  unknowns  are  updated  in  Eqs.  (13)  and  (15)  [23].  In  cases  where  no  prior 
structural  information  is  available,  a  pixel  basis  which  defines  a  set  of  regularly  spaced  pixels 
for  the  update  of  the  quantities  of  interest  is  used.  However,  when  spatial  priors  are  involved, 
it  is  crucial  to  ensure  that  enough  pixels  are  used  to  adequately  define  small  structural  regions. 
In  this  case,  a  set  of  regular  pixel  bases  are  introduced  in  each  region  of  interest,  which 
account  for  the  region’s  individual  shape  and  size.  A  semi-adaptive  method  which  allows  the 
number  of  pixels  in  each  region  to  be  selected  was  used  in  the  studies  presented  here.  The 
same  region-based  pixel  basis  is  used  throughout  the  iterative  reconstruction  process. 

3.  Methods 

Simulation  studies 

Test  domains  for  simulation  studies  were  derived  from  a  T1  weighted  MR  image  of  a  human 
breast  which  measured  approximately  10.5  cm  in  diameter.  A  2-D  slice  of  the  breast  volume 
was  discretized  into  a  mesh  of  approximately  2000  nodes  and  MR  image  intensity  thresholds 
were  used  to  assign  adipose  and  fibro-glandular  tissue  volumes  into  mesh  regions.  Figure  1 
shows  the  original  MR  breast  image  and  its  associated  discretized,  region-labeled  mesh.  The 
MR  image  originated  from  a  clinical  exam  with  our  MR  coupled  frequency  domain  NIR 
system  and  therefore  reveals  an  irregularly  shaped  breast  boundary  caused  by  the  fiber  optic 
array  slightly  compressing  the  breast  at  the  contact  positions.  A  cancerous  tumor  region  was 
added  as  a  target  anomaly  or  region  of  interest  for  these  studies  and  is  depicted  in  the  figure. 
A  second  test  case  with  a  tumor  region  located  near  the  center  of  the  domain  was  used  to 
explore  the  nonlinear  sensitivity  of  diffuse  optical  tomographic  techniques.  The  source- 
detector  positions  were  determined  directly  from  the  MR  image  and  represent  1 6  optical  fibers 
circumscribing  the  breast  domain  in  a  single  plane.  Light  is  detected  at  all  non-source  fiber 
positions  providing  a  total  of  240  measurements  of  intensity  and  phase  per  wavelength.  This 
simulated  configuration  matches  our  experimental  fluorescence  tomography  system. 


Fig.  1.  An  axial  T1  weighted  MR  slice  of  a  human  breast  is  shown  in  (a),  from  which  the  test 
domain  for  simulation  studies  was  derived.  Darker  regions  indicate  fibro -glandular  tissue 
imbedded  in  adipose  tissue  indicated  by  lighter  values.  The  image  was  acquired  during  a 
clinical  exam  with  one  of  our  MR-coupled  NIR  tomography  systems  and  shows  the 
indentations  caused  by  the  fiber  optic  probes.  The  test  domain  (b)  was  a  discretization  of  the 
MR  image  thresholded  into  regions.  The  yellow  anomaly  was  added  to  simulate  a  targeted 
cancerous  tumor. 

In  this  set  of  studies,  regions  were  assigned  values  in  terms  of  biologically  relevant  optical 
absorbing  (HbO,  dHb,  Water,  Lutetium  Texaphyrin)  as  well  as  scattering  parameters  (Mie 
scattering  amplitude  and  power).  These  values,  provided  in  Table  1,  represent  typical  known 
in  vivo  levels  of  endogenous  chromophores  and  are  consistent  with  previous  clinical  work 
[30,  31].  Though  it  is  unknown  exactly  how  Lutetium  Texaphyrin  (LuTex)  will  distribute  in  a 
human  breast,  the  concentrations  used  are  similar  to  those  reported  in  ex  vivo  studies  [17]  and 
at  least  represent  a  reasonably  complex  distribution  for  demonstrating  recovery  of 
fluorescence  yield.  Tissue  chromophore  concentrations  were  used  to  calculate  total  optical 
absorption  and  scattering  values  at  the  excitation  and  emission  wavelengths  based  on 
experimentally  determined  values  of  molar  extinction  coefficients.  Simulated  noisy  data  was 
generated  based  on  these  optical  properties  in  the  following  manner:  1)  frequency  domain 
data  for  a  light  source  at  the  excitation  wavelength,  2)  frequency  domain  data  for  a  light 
source  at  the  emission  wavelength,  and,  3)  Continuous  Wave  (CW)  fluorescence  emission 
data  at  the  emission  wavelength.  This  represents  a  total  of  3  data  sets  for  a  given  imaging 
session:  two  frequency  domain  data  sets  for  determining  background  optical  properties  and 
one  CW  fluorescence  emission  data  set.  CW  fluorescence  emission  data  was  used  to  match 
the  capabilities  of  our  experimental  system  and  results  in  a  simplification  of  Eq.  (2)  which  can 
be  handled  by  setting  co  =  0. 


Table  1.  Chromophore  concentrations  and  scattering  parameter  values  assigned  to  the  mesh 
regions  in  the  simulation  studies 


Oxy¬ 

hemoglobin 

(MM) 

Deoxy¬ 

hemoglobin 

(mM) 

Water 

(%) 

Lutetium 

Texaphyrin 

(HM) 

Scattering 

Amplitude 

Scattering 

Power 

Adipose 

10 

10 

50 

0.3 

1.0 

1.0 

Fibro- 

15 

15 

60 

0.5 

1.1 

1.1 

glandular 

Tumor 

20 

18 

90 

1.0 

1.2 

1.15 

Given  the  modest  Stoke ’s  shift  of  LuTex,  depicted  in  Fig.  2,  the  choice  of  excitation 
wavelength  has  practical  implications  for  experimental  work.  The  difficulty  in  filtering  the 


excitation  signal  precludes  the  use  of  Lutex’s  NIR  absorption  peak  (about  735  nm)  to  excite 
the  fluorophore.  In  this  study,  a  690  nm  excitation  wavelength  was  used  for  both  simulated 
and  experimental  data  acquisition.  In  addition  to  normally  distributed  random  noise  added  to 
the  frequency  domain  data  (5%  amplitude  and  1  degree  phase)  and  CW  fluorescence  emission 
(10%  intensity),  excitation  signal  leakage  through  the  filter  was  added  to  the  CW  fluorescence 
emission  intensity  to  simulate  a  typical  7  OD  rejection  of  the  excitation  intensity.  This 
number  comes  directly  from  experimentally  measured  rejection  estimates  for  the  filters  in  the 
tomography  system  used  in  the  experiments  performed  here. 


Fig.  2.  The  normalized  absorbance  and  fluorescence  emission  spectra  of  Lutetium  Texaphyrin 
are  shown. 

The  general  image  reconstruction  protocol  was  as  follows: 

1)  Reconstruct  for  optical  properties  at  the  excitation  wavelength,  jaax  and  jlisx’,  with 
frequency  domain  data, 

2)  Reconstruct  for  optical  properties  at  the  emission  wavelength,  juam  and  jusm’,  with 
frequency  domain  data  collected  using  a  laser  source  at  the  emission  wavelength, 

3)  Use  the  reconstructed  optical  properties  and  fluorescence  intensity  data  to  recover 
fluorescence  yield. 

The  same  reconstruction  algorithm  was  used  to  determine  background  optical  properties  in 
steps  (1)  and  (2)  and  is  based  on  previously  reported  work  [32,  33].  Initial  estimates  for  all 
parameters  were  generated  using  homogenous  fitting  algorithms  which  enforce  a  single  value 
for  all  nodes.  The  Jacobian  matrix  was  calculated  on  a  fine  mesh  of  approximately  2000 
nodes  and  interpolated  onto  a  course  reconstruction  pixel  basis  for  inversion.  A  30  by  30 
pixel  reconstruction  basis  was  used  for  the  no-priors  case.  The  spatial  priors  reconstruction 
used  a  newly  developed  semi-adaptive  pixel  basis  that  redistributes  the  pixels  based  on  the 
region  information,  as  described  in  section  2.2.3.  This  method  ensures  that  each  region 
contains  an  adequate  number  of  nodes  to  approximate  the  internal  structure  of  the  domain. 
Convergence  was  defined  as  less  than  a  2%  change  in  projection  error  between  successive 
iterations  for  the  frequency  domain  optical  properties  algorithm  and  less  than  a  1%  change  in 
projection  error  between  iterations  for  the  fluorescence  yield  reconstruction.  Similar 
algorithmic  parameters  were  used  for  the  phantom  studies  and  are  described  in  further  detail 
below. 

Phantom  studies 

A  spectrometer  based  tomographic  imaging  system  which  couples  directly  into  a  Philips  3T 
MRI  magnet  was  used  to  acquire  fluorescence  emission  spectra.  The  system,  depicted  in  Fig. 


3,  is  composed  of  16  spectrometers  with  low  noise  CCD  cameras  cooled  to  -70  C.  Sixteen  13 
meter  long  silica  fiber  bundles  circumscribe  the  imaging  domain  in  contact  mode  and  couple 
into  the  spectrometers  via  custom  designed  input  optics  with  a  six  position  automated  filter 
wheel.  Fluorescence  emission  signals  are  processed  with  long  pass  interference  filters  before 
entering  the  spectrometer.  Each  fiber  bundle  contains  eight  400  jam  fibers,  seven  of  which 
directly  couple  the  tissue  surface  to  the  spectrometer  input  while  the  eighth  couples  the  tissue 
surface  to  the  light  source.  This  detector  configuration  minimizes  coupling  losses  and 
provides  parallel  detection  of  full  spectra  for  each  source  position.  Inter-fiber  tissue  coupling 
calibration  factors  were  determined  prior  to  imaging  using  a  cylindrical  homogenous  phantom 
with  a  centrally  located  source/detector.  The  custom  manufactured  input  optics  were  designed 
to  optimize  filtering  and  maximize  the  light  detection  efficiency  of  the  spectrometers.  Camera 
exposure  times  can  be  adjusted  to  the  detected  light  level  to  maximize  the  signal  to  noise  ratio. 
With  15  detectors  per  16  source  positions,  a  total  of  240  measurements  are  recorded  for  a 
given  acquisition. 


Fig.  3.  The  experimental  spectrometer-based  system  depicted  at  left  couples  directly  into  the 
MR  via  13  meter  fiber  optic  bundles.  Sixteen  spectrometers  are  computer  controlled  for  rapid 
image  acquisition  (left  photograph).  An  animal  interface  (right  photograph)  is  composed  of  a 
rodent  MR  coil  custom  built  by  Philips  Research  Hamburg  to  accommodate  the  optical  fiber 
array  for  simultaneous  MR  and  NIR  fluorescence  imaging. 

Lutetium  Texaphyrin  provided  by  Pharamcyclics  was  diluted  in  water  and  used  as  the 
imaging  fluorophore.  The  test  domain  was  a  solid  5.5  cm  diameter  hardened  epoxy  phantom 
with  scatterer  and  absorbers  created  by  titanium  dioxide  powder  and  India  ink,  respectively 
[34,  35].  The  phantom  had  a  14  mm  hole  located  approximately  12  mm  from  the  phantom 
center.  The  background  optical  properties  were  jnax  =  0.005  and  jnsx’  =  1.0  mm"1,  measured 
with  a  frequency  domain  system  near  the  excitation  wavelength.  Unlike  the  simulation 
experiments,  the  optical  properties  were  assumed  constant  throughout  the  domain  in  this 
experiment.  These  values  were  also  used  as  the  optical  properties  at  the  emission  wavelength. 
The  hole  was  filled  with  a  solution  of  1%  Intralipid  to  match  the  scattering  value  of  the 
background,  and  varying  concentrations  of  Lutex  (0.3125  jliM  to  5  jliM)  were  added.  This 
represents  a  simple  test  case  for  investigating  the  imaging  response  to  varying  concentrations 
of  fluorophore.  The  excitation  source  was  a  690  nm  laser  diode  which  matches  the 
wavelength  used  in  the  simulation  studies.  Total  acquisition  time  for  the  fluorescence 
emission  was  less  than  4  minutes  (total  of  240  data  points). 

Even  after  processing  the  collected  light  with  a  720  nm  long  pass  interference  filter 
(Omega  Optics)  which  provides  7  OD  rejection  of  the  excitation  light  as  well  as  the  filtering 
offered  by  the  spectrograph  grating,  emission  spectra  recorded  by  the  detector  are  composed 
of  a  sum  of  the  pure  fluorescence  signal  and  excitation  bleed-through.  In  order  to  decouple 
these  signals,  previously  recorded  “basis  spectra”  of  the  bleed-through  signal  and  the  pure 
fluorescence  signal  are  fit  to  the  data.  The  process  is  illustrated  in  Fig.  4  for  one  measurement 


at  a  single  detector.  A  linear  least  squares  algorithm  operating  on  the  basis  spectra  quantifies 
the  amount  of  excitation  bleed-through  and  true  fluorescence  signal  that  exists  in  each  spectral 
recording,  further  reducing  the  effect  of  excitation  bleed-through.  This  is  accomplished  by 
minimizing  the  summation 

S  =  f}y,-{aF^i)  +  bG(Ai))]2  (16) 

i= 1 

with  respect  to  a  and  b ,  where  yt  is  the  measured  intensity  at  a  given  wavelength  pixel,  F  and 
G  are  the  excitation  and  fluorescence  basis  spectra,  a  and  b  are  the  coefficients  recovered  in 
the  minimization  procedure,  and  N  is  the  number  of  wavelength  pixels  per  spectrum.  Once 
fit,  the  fluorescence  signal  is  integrated  and  becomes  the  fluorescence  emission  intensity  data 
for  the  reconstruction  algorithm. 


(a)  (b) 


CCD  Pixel  number  CCD  Pixel  number 


Fig.  4.  An  example  of  a  pair  of  basis  spectra  for  the  excitation  and  fluorescence  light  (in 
counts/s  as  a  function  of  CCD  pixel  number)  is  shown  in  (b).  These  spectra  are  recorded  for 
each  detector  prior  to  imaging.  In  practice,  the  basis  spectra  are  used  to  perform  a  least  squares 
fit  (a)  to  the  spectrum  measured  for  each  source-detector  pair  to  determine  the  relative 
contribution  of  the  fluorescence  and  excitation  light  to  the  measured  response. 

To  investigate  the  improved  imaging  capability  of  reconstructing  data  with  spatial  priors, 
each  data  set  was  reconstructed  with  and  without  spatial  “soft”  priors.  Since  the  domain  was 
easily  characterized  geometrically,  spatial  prior  information  was  determined  by  direct  manual 
measurement  of  the  phantom.  This  information  was  encoded  in  the  fine  mesh  used  for 
reconstruction.  Since  background  optical  properties  were  previously  determined  with  our 
frequency  domain  system,  they  were  assumed  known  for  image  reconstruction.  Convergence 
was  defined  as  less  than  a  1  %  change  in  projection  error  between  iterations.  All  images  were 
reconstructed  using  a  2GHz  Centrino  Duo  laptop  with  2GB  RAM  running  Windows  XP. 

4.  Results 

Simulation  results 

Figures  5  and  6  display  the  target  values  of  the  two  test  cases  along  with  images  reconstructed 
using  no  priors  and  spatial  soft  prior  information.  Qualitatively,  image  recovery  using  spatial 
priors  produces  significantly  more  accurate  images.  Spatial  priors  preserve  the  general 
internal  structure  of  the  Lutex  distribution,  detail  that  is  lost  almost  entirely  in  the  no-priors 
case.  Images  of  fluorescence  yield  show  the  most  dramatic  difference  between  the  spatial 
priors  and  no  priors  reconstructions.  Without  spatial  priors,  the  algorithm  appears  to  have  no 
ability  to  recover  “cancer”  regions  of  elevated  fluorescence  yield  for  these  complicated  cases. 
However,  incorporating  spatial  priors  results  in  images  that  qualitatively  appear  accurate  and 
quantitatively  are  reasonably  close  to  the  true  values.  Figure  7  provides  1-D  cross-sections 
near  the  y-axis  of  each  domain.  These  plots  confirm  an  inability  of  the  no-priors  imaging 
algorithm  to  recover  the  simulated  fluorescent  tumor  in  either  test  field.  Alternatively,  the 


spatial  prior-based  imaging  algorithm  not  only  picks  out  the  objects  of  interest,  but  provides 
fairly  accurate  reproductions  of  the  complicated  structure  of  simulated  fibro-glandular  layers. 
Mean  values  for  the  simulated  cancer  region  near  the  domain  center  are  79%  and  51%  of  the 
true  values  for  the  spatial  priors  and  no-priors  reconstructions,  respectively.  These  numbers 
change  to  75%  and  45%  for  the  case  with  the  cancer  region  closer  to  the  edge  of  the  domain. 
They  represent  a  significant  improvement  in  imaging  performance;  however,  they  alone  do 
not  illustrate  the  full  impact  of  incorporating  spatial  priors.  The  cross-sectional  plots  indicate 
that  the  no-priors  images  contain  virtually  no  spatial  discrimination  of  the  cancer  regions. 
Regional  contrasts  are  depicted  in  Table  2  and  further  illustrate  a  dramatic  overall 
improvement  in  cancer  region  quantification  with  spatial  priors. 


Target  distribution  to  No  spatial  priors  Spatial  “soft”  priors 

generate  data  reconstructed  image  reconstructed  image 


Fig.  5.  Target  and  recovered  values  of  pa,x,  ps,x,  pa,m,  p  s,m,  and  fluorescence  yield,  ri|Liaf ,  for 
reconstruction  implementations  using  no  prior  information  and  with  spatial  prior  information. 
In  this  case,  the  simulated  cancer  region  is  near  the  edge,  which  is  known  to  be  easier  to 
recover  without  spatial  priors.  The  image  scales  are  at  right. 


Fig.  6.  Target  and  recovered  values  of  pa,x,  |l4,x,  pa,m,  p^m,  and  fluorescence  yield,  ri|Liaf ,  for 
reconstruction  implementations  using  no  priors  and  spatial  priors.  In  this  case,  the  simulated 
tumor  region  is  near  the  center  of  the  imaging  domain,  which  is  known  to  be  more  difficult  to 
recover  accurately.  Image  scales  are  at  right. 


Table  2.  Target  and  recovered  fluorescence  yield  regional  contrasts  for  the  images  in  Figs.  4 
and  5. 


Shallow  tumor 

Target  Contrast 

No  spatial  priors 

Spatial  “soft”  priors 

Tumor  :  fibro- 

2.0  :  1 

1.1  :  1 

1.6  :  1 

glandular 

Tumor  :  adipose 

3.3  :  1 

1.3  :  1 

2.3  :  1 

Deep  tumor 

Target  Contrast 

No  spatial  priors 

Spatial  “soft”  priors 

Tumor  :  fibro- 

2.0  :  1 

1.2  :  1 

1.8  :  1 

glandular 

Tumor  :  adipose 

3.3  :  1 

1.4  :  1 

2.4  :  1 

Reported  ratios  were  calculated  using  the  mean  values  in  each  region. 


Fig.  7.  Cross  sectional  plots  of  fluorescence  yield  are  shown  for  the  simulated  imaging 
domains  in  (a)  the  case  with  an  object  near  the  edge  and  (b)  the  case  with  an  object  near  the 
center.  In  both  cases,  the  cross  section  is  in  the  y-direction  just  off  center  from  x  =  0.  The 
solid  line  represents  the  target  value,  the  small  dotted  line  the  recovered  value  using  a  no-priors 
based  algorithm,  and  the  dashed  line  the  recovered  value  using  spatially  guided  reconstruction. 


In  addition  to  improved  qualitative  and  quantitative  accuracy,  the  spatial  prior  algorithm 
reduced  the  reconstruction  time  significantly.  The  full  reconstruction  time  for  the  no-priors 
case,  including  background  optical  property  estimation,  was  just  under  9  minutes  for  both  the 
central  and  superficial  tumor  cases.  These  times  were  reduced  to  less  than  3  minutes  and  90 
seconds  using  spatial  priors.  In  both  cases,  structural  information  guided  the  algorithm  to  a 
convergent  solution  in  far  fewer  iterations  than  the  no-priors  cases. 

Phantom  results 

Figure  8  shows  fluorescence  yield  images  recovered  from  phantom  data  using  both  no-prior 
and  spatial  prior  image  reconstruction  approaches.  A  qualitative  assessment  of  the  images 
reveals  a  dramatic  benefit  of  the  spatial  prior  on  image  formation.  The  fluorescent  object’s 
borders  are  more  clearly  defined  for  all  concentrations  of  Lutex  when  using  spatial  priors. 
Furthermore,  the  values  of  fluorescence  yield  throughout  the  region  of  interest  are  more 
homogeneous,  and  therefore,  more  similar  to  the  actual  distribution  for  spatially  guided 
reconstructions.  Incorporating  spatial  priors  also  suppresses  edge  artifacts  significantly.  This 
is  most  apparent  in  images  of  phantoms  with  low  Lutex  concentration.  Figure  9  shows  a  full 
scale  image  of  the  0.3125  pM  phantom.  Artifacts  near  the  boundary  virtually  dominate  the 
no-priors  image  and  represent  the  largest  fluorescent  yield  values  in  that  imaging  domain. 
None  of  these  artifacts  exist  in  the  spatially  guided  image  generated  from  the  same 
fluorescence  emission  data.  The  spatially  guided  image  also  includes  an  easily  discernable 


fluorescent  object  at  the  correct  location  which  is  difficult  to  define  in  the  no-priors  case. 
These  results  indicate  that  for  low  fluorophore  concentrations  in  this  imaging  domain,  the 
correct  fluorescent  object  would  be  identified  only  if  prior  structural  information  were 
incorporated  in  the  reconstruction. 


Fig.  8.  Recovered  images  of  fluorescence  yield  are  shown  for  varying  concentrations  of  Lutex. 
The  14  mm  diameter  fluorescent  inclusion  was  embedded  in  a  55  mm  diameter  solid  epoxy 
tissue  simulating  phantom.  Images  were  generated  from  the  same  data  using  algorithms  based 
on  no-priors  and  spatial  soft  prior  implementations. 


Fig.  9.  A  narrower  colorbar-scale  version  of  the  0.3125  jaM  Lutex  phantom  images  shown  in 
Fig.  8  further  illustrates  the  improvement  in  image  accuracy  for  the  spatially  guided  algorithm. 

Both  techniques  suffer  from  a  bias  in  the  recovered  fluorescence  yield  which  results  in  a 
positive  value  in  the  region  of  interest  for  the  case  with  no  fluorophore.  This  is  most  likely 
caused  by  bleed-through  of  the  excitation  light  in  the  emission  measurements  coupled  with  an 
imperfect  spectral  fitting  technique  which  results  in  positive  values  for  fluorescence  emission 
intensity  even  in  phantoms  containing  no  Lutex  fluorophore.  Incorporating  additional 
filtering,  such  as  band-pass  source  filtering,  and  optimizing  the  spectral  fitting  routine  should 
address  the  bias  signal  and  decrease  the  potential  for  incorrect  quantification. 

In  addition  to  improving  image  quality  and  fluorescence  yield  quantification,  spatial  prior- 
based  reconstructions  converged  in  substantially  less  time  than  the  no-priors  algorithms. 
Reconstruction  times  ranged  from  50  -  90  seconds  when  not  using  spatial  priors  and 
approximately  22  seconds  for  spatially  guided  reconstructions,  marking  an  improvement  of 
just  over  75%  in  some  cases.  These  numbers  represent  only  the  fluorescence  yield 
reconstructions  and  not  the  recovery  of  background  optical  properties  since  jua  and  jlis’  were 
assumed  known  from  prior  measurements.  In  cases  requiring  the  recovery  of  background 
optical  properties,  the  improvement  in  reconstruction  time  will  be  similar  to  that  previously 
quoted  in  the  simulation  results. 

5.  Discussion 

This  study  introduced  an  effective  method  to  incorporate  MR-derived  tissue  morphology  for 
imaging  fluorescence  yield  at  depth  in  tissue.  The  conceptual  assertion  is  that  diffuse 
tomography  will  be  more  successful  as  an  imaging  modality  when  combined  with  pre-existing 
imaging  systems  which  have  higher  spatial  resolution,  such  as  MRI.  Simulation  and  phantom 
studies  were  used  here  to  validate  the  method  prior  to  implementation  as  a  full  scale  system. 
Using  Lutex  phantom  data,  the  soft  spatial  prior  implementation  improved  qualitative  and 
quantitative  imaging  response  of  fluorescence  yield.  Prior  information  also  suppressed  image 
artifacts  and  more  accurately  represented  the  internal  distribution  of  fluorophore.  In  the  no¬ 
priors  case,  edge  artifacts  dominate  the  image  for  lower  concentrations  of  fluorophore,  and 
provide  a  misleading  interpretation  of  the  internal  distribution  of  the  fluorescent  agent.  The 
phantom  studies  represented  simple  distributions  of  optical  properties  and  fluorophore 
concentration  and  improvements  in  image  quality  were  still  substantial.  It  is  expected  that 
soft  prior  implementations  will  benefit  imaging  performance  to  an  even  greater  extent  in 
complex  tissue  domains,  an  assertion  bom  out  in  the  simulation  studies. 

A  simulated  breast  mesh  derived  directly  from  an  MR  image  of  a  human  breast  served  as 
the  test  bed  for  a  complex  tissue  domain.  The  simulations  demonstrated  no  ability  to  recover 
the  internal  distribution  of  the  complicated  domain  without  the  use  of  spatially  guided 
reconstmctions.  Qualitatively,  fluorescence  yield  images  generated  without  spatial  priors  had 
little  resemblance  to  the  target  domain  and  completely  disregarded  the  18mm  simulated 
cancer  region,  as  evidenced  by  1  -D  cross-sectional  plots  of  fluorescence  yield.  Breakdown  of 
the  images  in  the  no-priors  case  is  likely  due  to  poorly  recovered  background  optical 
properties  as  well  as  the  complexity  of  the  fluorescence  yield  distribution  itself.  Certainly, 


part  of  the  improvement  in  fluorescence  yield  image  accuracy  can  be  attributed  to  improved 
images  of  background  optical  properties. 

In  these  studies,  it  was  assumed  that  the  fluorophore  distribution  correlates  directly  to  the 
fatty,  fibro-glandular,  and  tumor  tissue  layers,  which  themselves  exhibit  no  significant  intra¬ 
region  heterogeneity,  in  a  given  imaging  domain.  While  this  assumption  appears  reasonable 
based  upon  the  biology  of  the  tissue  for  endogenous  chromophores  [30,  31],  further  studies  of 
the  uptake  of  Lutex  in-vivo  are  required  to  determine  a  correlation  between  tissue  type  and 
fluorophore  localization.  Should  this  assumption  prove  to  be  unrealistic,  the  “soft”  spatial 
prior  approach  offers  some  latitude  in  terms  of  correctly  identifying  structural  prior 
information.  The  algorithmic  implementation  groups  tissue  regions  together  in  a  region- 
specific  regularization  and  allows  individual  nodes  in  those  regions  to  update  independently. 
As  opposed  to  “hard”  prior  approaches  where  nodal  values  in  a  given  region  are  assumed 
homogenous,  a  soft  prior  technique  may  recover  positive  fluorescent  objects  not  directly 
encoded  in  the  spatial  prior  information.  This  is  a  subject  of  further  investigation. 

The  experimental  system  introduced  here  couples  directly  into  a  Philips  3T  MRI  to 
provide  simultaneous  MR  and  NIR  fluorescence  imaging.  Simultaneous  imaging  simplifies 
co-registration  of  the  MR  image  with  the  NIR  reconstruction  domains  and  reduces  overall 
acquisition  time.  Optimization  of  this  system  for  in  vivo  imaging  is  underway  for  both  small 
animal  and  human  breast  imaging. 

6.  Conclusion 

A  spatially  guided  image  reconstruction  implementation  based  on  prior  knowledge  of  tissue 
morphology  was  shown  to  provide  significant  improvements  in  fluorescence  yield  recovery  in 
complicated  tissue  volumes  and  to  be  highly  beneficial  for  simple  domains.  Specifically,  both 
phantom  and  simulation  results  demonstrated  dramatic  improvements  in  recovery  and 
quantification  of  features  in  the  fluorescence  distribution.  Structural  guidance  also  reduced 
image  reconstruction  time  substantially.  A  newly  developed  experimental  system  couples 
full-volume  deep-tissue  fluorescence  spectroscopy  capabilities  directly  into  the  MR  bore  for 
simultaneous  MR-NIR  fluorescence  image  acquisition.  The  results  presented  here  show 
promise  for  this  approach  in  all  cases  and  tissue  volumes  considered.  Extensions  of  this  study 
are  underway  to  determine  the  effect  of  incorrectly  identifying  the  structural  prior,  especially 
in  cases  where  the  MR  images  produce  false  negative  or  false  positive  readings.  A  full 
characterization  of  the  imaging  limits  of  the  experimental  system  will  further  complement  the 
results  presented  here.  It  is  clear  that  incorporating  anatomical  features  derived  from  MR 
images  in  DOFT  image  reconstruction  will  improve  sensitivity  to  lower  concentrations  of 
fluorophore,  qualitative  accuracy,  and  fluorescence  yield  quantification  in-vivo. 
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